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Abstract 

Parallel MRI (pMRI) and compressed sensing MRI (CS-MRI) have been considered as two distinct reconstruction 
problems. Inspired by recent k-space interpolation methods, an annihilating filter based low-rank Hankel matrix 
approach (ALOHA) is proposed as a general framework for sparsity-driven k-space interpolation method which 
unifies pMRI and CS-MRI. Specifically, our framework is based on the fundamental duality between the transform 
domain sparsity in the primary space and the low-rankness of weighted Hankel matrix in the reciprocal space, 
which converts pMRI and CS-MRI to a k-space interpolation problem using structured matrix completion. Using 
theoretical results from the latest compressed sensing literatures, we showed that the required sampling rates for 
ALOHA may achieve the optimal rate. Experimental results with in vivo data for single/multi-coil imaging as well 
as dynamic imaging confirmed that the proposed method outperforms the state-of-the-art pMRI and CS-MRI. 
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I. Introduction 


Magnetic resonance imaging (MRI) is an imaging system that sequentially acquires k-space data corresponding to 
the Fourier transform of an object. This enables us to apply various advanced signal processing techniques. Recently, 
compressed sensing theory Q, Q has been used in accelerated MRI Q-Q. Compressed sensing algorithms can 
restore original signals from much less k-space data by exploiting the sparsity of an unknown image in total variation 
(TV) or wavelet transform domains, and incoherent sampling schemes such as Gaussian random or Poisson disc 
are usually required. Accurate MRI reconstruction from less data makes compressed sensing a hot topic in the 
research community; thus, it has been applied across many different application areas such as in pediatric imaging 
Q, dynamic cardiac MRI Q-Q, perfusion imaging 110|, angiography pT| , and so on. 

On the other hand, parallel MRI (pMRI) |[T^-| 141 exploits the diversity in the receiver coil sensitivity maps 
that are multiplied by an unknown image. This provides additional spatial information for the unknown image, 
resulting in accelerated MR data acquisition through k-space sample reduction. Representative parallel imaging 


algorithms such as SENSE (sensitivity encoding) |12| or GRAPPA (generalized autocalibrating partially parallel 


acquisitions) require regularly sampled k-space data for computationally efficient reconstruction. Moreover, 
additional k-space data, the so-called auto calibration (ACS) lines, are often required to estimate the coil sensitive 


maps or GRAPPA kernels |T^. 

Because the aim of the two approaches is accelerated acquisition by reducing the k-space data, extensive research 
efforts have been made to synergistically combine the two for further acceleration. One of the most simplest 
approaches can be a SENSE type approach that explicitly utilizes the estimated coil maps to obtain an augmented 
compressed sensing problem: 


mm ||VFf 111 


subject to g = 


gl 


■A[Ni]- 

.S'-. 




( 1 ) 


where f and gj denote the unknown image and the k-space measurements from the i-th coil, respectively; A is 
a subsampled Eourier matrix; PF is a sparsifying transform, and [5*] denotes a diagonal matrix whose diagonal 
elements come from the f-th coil sensitivity map. The multichannel version of k-t EOCUSS 0 is one of the 
typical examples of such approaches. On the other hand, (i-SPIRiT (fi- iTerative Self-consistent Parallel Imaging 


Reconstruction) 1151 utilizes the GRAPPA type constraint as an additional constraint for a compressed sensing 
problem: 


mm 

F 


I^'FI 


1,2 


subject to G = AF 

Vec(F) = M ■ Vec(F) 

where || • ||i ^2 denotes the (l,2)-mixed norm of a matrix, F = [fi £2 • • • f^] , G = [gi g 2 


( 2 ) 

(3) 

(4) 
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the unknown images and their k-space measurements for the given set of coils, and tk denote a discrete wavelet 
transform matrix, and M is an image domain GRAPPA operator, and Vec(-) is the vectorization operator. In both 
approaches, an accurate estimation of coil sensitivity maps or GRAPPA kernel is essential to fully exploit the coil 
sensitivity diversity. 

In order to overcome these difficulties, calihration-less parallel imaging methods have been extensively investi¬ 


gated, among which SAKE (simultaneous autocalibrating and k-space estimation) 1161 represents one of the first 
steps. In SAKE, the missing k-space elements are reconstructed by imposing the data consistency and the structural 
maintenance constraints of the block Hankel structure matrix. However, the origin of the low rankness in the Hankel 
structured matrix for the case of a single coil measurement was not extensively investigated, and it was not clear 
whether SAKE could outperform the compressed sensing approach when it is applied to single coil data. Haidar 
|T7), |T8) later discovered that a Hankel structured matrix constructed by a single coil k-space measurement is 
low-ranked when an unknown image has finite support or a slow-varying phase. Based on this observation, he 
developed the so called LORAKS (Low-rank modeling of local k-space neighborhoods) algorithm p7| and its 
parallel imaging version, P-LORAKS (Low-rank modelling of local k-space neighborhoods with parallel imaging 


data) 118|. However, it is not clear how the existing theory can deal with large classes of image models that are 
not finite supported but can be sparsified using various transforms such as wavelet transforms or total variations 
(TV), etc. 

Therefore, one of the main goals of this paper is to develop a theory that unifies and generalizes k-space 
low-rank mefhods fo also allow for fransform sparsify models which are critical for pracfical MRI applicafions. 
Toward fhis goal, we show fhaf fhe fransform domain sparsity in the signal space can be directly related to the 
existence of annihilating filters m in the weighted k-space. Interestingly, the commutative relation between 
an annihilating filter and weighted k-space measurements provides a rank-deficient Hankel structured matrix, whose 
rank is determined by the sparsity level of the underlying signal in a transform domain. Therefore, by performing 
a low-rank matrix completion approach, the missing weighted k-space data in the Hankel structured matrix can be 
recovered, after which the original k-space data can be recovered by removing the weights. 

Interestingly, the new framework is so general that it can cover important compressed sensing MR approaches in 
very unique ways. For example, when an image can be sparsified wifh wavelef fransforms, fhe low rank sfrucfured 
matrix completion problem can be solved using a pyramidal decomposition after applying scale dependent k-space 
weightings. Specifically, a low rank Hankel mafrix completion algorifhm can be progressively applied from fine fo 
coarse scale to reduce the overall computational burden while maintaining noise robustness. In addition, we show 
that there exist additional inter-coil annihilating filter relationships that are unique in pMRI, which can be utilized 
to construct a concatenated Hankel matrix which is low ranked. We further substantiate that the multi-channel 
stacking of the weighted Hankel structure matrix may fully exploit the coil diversity thanks to the relationship to 
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the algebraic bound of multiple measurement vector (MMV) compressed sensing \22\-\24\. 


Another important advantage of the proposed algorithm is that, compared to the existing CS-MRI, the recon¬ 
struction errors are usually scattered throughout the entire images rather than exhibiting systematic distortion along 
edges because the annihilating filter relationships are specifically designed for estimating the edge signals. Given that 
many diagnostic errors are caused by the systematic distortion of images, we believe that our ALOHA framework 
may have great potential in clinical applications. 

The remainder of this paper is as follows. Section discusses the fundamental duality between the transform 
domain sparsity and the low-rankness of Hankel structured matrix in k-space. In Section a detailed description 
of the proposed method for MR reconstruction will be provided. Section then explains the implementation detail. 
Experimental results are provided in Section [Vj which is followed by the discussion in Section VI and conclusion 
in Section IVlIl 

II. Fundamental Duality between Sparsity and Low-Rankness 

This section describes the fundamental dual relationship between transform domain sparsity and low rankness 
of weighted Hankel matrix, which is the key idea of the proposed algorithm. For better readability, the theory here 
is outlined by assuming 1-D signals, but the principle can be extended for multidimensional signals. 

Note that typical signals of our interest may not be sparse in the image domain, but can be sparsified in a 
transform domain. For example, consider L-spline signal model from the theory of sparse stochastic processes 
|. Specifically, fhe signal / of our inferesf is assumed to satisfy the following partial differential equation: 


Lf = w 


(5) 


where L denotes a constant coefficient linear differential equation (or whitening operator in |25|, |26|): 

L := axd^ -f aK-id^~^ -|- ... + aid + oq (6) 


and ru is a driving continuous domain sparse signal given by 

fc-i 

j=0 


(7) 


Here, without loss of generality, we set r = ni for a positive integer ni G Z. This model includes many class of 


signals with finite rate of innovations 119|-121|. For example, if the underlying signal is piecewise constant, we 
can set L as the first differentiation. In this case, / corresponds to the total variation (TV) signal model, and this 
TV signal model will be extensively used throughout the paper. 

Now, by taking the Fourier transform of Q, we have 


k-l 


y{uj) := J='{Lf{x)} = i{uj)f{uj) = '^aje~ 

j=0 


( 8 ) 
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where 


l{uj) = axiiujr + ^ + ... + ai{iuj) + oq 


(9) 


In standard Nyquist sampling, we should measure discrete set of Fourier samples at {wj}™ ^ from a deterministic 
grid, whose grid size should he set to the Nyquist limit A = 27r/ni to avoid aliasing artifacts; so the discrete 
specturm can he represented as 


fc-i 


( 10 ) 


y[m] := = [[mA]/[mA] = 

j=0 

for m G [0, • • • ,ni — 1]. The discrete spectral sampling model in Eq. ( [T0| ) implies that the unknown signal in the 
image domain is an infinite periodic streams of Diracs with a period rei, which is indeed a signal with the finite rate 


of innovation (FRI) with rate p = 2k/ni \ 19|-|211. Therefore, theoretical results from the FRI sampling theory can 
he used p^-|[2H. In particular, the FRI sampling theory tells us that we can find a minimum length annihilating 
filter h[n] such that 


{h*y)[n] = '^h[l]y[n-l] = Q, Vn. 


( 11 ) 


1=0 


The specific form of fhe minimum length annihilafing filter h[n] for the case of ( fTO] ) can he found in 119|, which 
has the following z-transform 


k-l 


h{z) = 




( 12 ) 


1=0 


j=0 


whose filter length is given hy k + 1 119|. 


Now, one of the most important advances of the proposed method over the existing sampling theory of FRI 
signals is that, rather than using the minimum length annihilating filter, we allow a longer length annihilating filter 
which is essential to make the associated Hankel matrix low-ranked. Specifically, if h[n] is a minimum length 
annihilating filter with k + 1 filter taps, then for any ki>l tap filter a[n], it is easy to see that the following filter 
with K = k + ki taps is also an annihilating filter for y[n]\ 


ha[n] = (d * h)[n] . 

Accordingly, the matrix representation of (ha * y) [n] = 0 is given hy 

‘^(y)ha = 0 


(13) 


where ho denotes a vector that reverses the order of the elements in 


ha = 


^a[0], • • • ,/la[K - 1] 


1 T 


(14) 
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and 


'^(y) 


y[-i] 

y[o] 

y[K - 2] 

y[o] 

m 

y[K- 1] 

y[i] 

y[2] 

y[>^] 

y[ni - k] 

y[ni-fi: + l] ••• 

y[n - 1] 

y[ni - K + 1] 

y[ni - k + 2] ■■■ 

y[ni] 





(15) 


Accordingly, by choosing ni such that ni — k + 1 > k and removing the boundary data outside of the sample 
indices [0, • • • , ni — 1], we can construct the following matrix equation: 


Jf'(y)ha = 0, 


(16) 


where the Hankel structure matrix J^{y) G k+i)xk constructed as 


^{y) 


y[o] y[i] 


y[ni - k] y[ni - k + 1] 


Then, we have the following result: 


y[K - 1] 


y[ni - 1] 


(17) 


Proposition 2.1. Let k denotes the number of Diracs within the support [0,ni]. Suppose, furthermore, the annihi¬ 
lating filter length n is given by ni — k 1 > k. Then, for a given Hankel structured matrix J^{y) in ( |17[ ), we 
have 


RANKjf'(y) < k, 


(18) 


where RANk(-) denotes a matrix rank. 


Proof Note that can be represented as 


ha = ‘^(h)a 


(19) 


where a = [d[0] • • • d[ki — 1]] and ^(h) G is a Toeplitz structured convolution matrix from h: 



r h[ 0 ] 

0 

0 


h[l] h[0] 

0 

^(h) = 

h[k — 1] h[k — 2] 

• • • h[k — ki] 


0 

0 

h[k — 1] _ 


where k = k ->r ki. Since '^(h) is a convolution matrix, it is full ranked and we can show that 


( 20 ) 


dim'^(h) = ki 













7 


where dim(-) denotes the dimension of a matrix. Moreover, the range space of '^(h) now generates the null space 
of the Hankel matrix, so it is easy to show 

ki = dim^(h) < dimNULjf'(y), 

where Nul(-) represent a null space of a matrix. Hence, we have 

RANK,^(y) = min{K,ni — ki + 1} — dimNUL,^(y) < k — ki = k. (21) 

Q.E.D. □ 

This rank condition can he tightened under more restricted set-ups. Specifically, ifm — k + 1 > k + 1, the authors 
in | |271 l derived the following decomposition: 

J^{y) = L D , 


( 22 ) 


where 


L = 


R = 


1 

1 

1 




zq 

Zi 

• • • Zk- 

-1 


(23) 

ni-K 

0 

ni-K 

■■■ 

— K 

1 . 



1 

1 

1 ■ 




ZQ 

Zl 

• • Zk-1 

E 

(24) 

K-1 

0 


■ ■ z'^~^ 





and 


D = 


Co 0 

0 Cl 


0 

0 


(25) 


0 0 ••• Cfc_i_ 

and Zj = for j = 0, • • • , A: — 1. Because L and R are Vandermonde matrices that are full column-ranked, 

in this case we can further show that 

RANK^(y) = RANK(T») = k. 

However, as shown in | [28| , Proposition |2. 1 1 covers more general classes of Hankel matrices that do not necessarily 


have the Vandermonde decomposition structure in (22). For example. Proposition |2.1| is still valid even when the 


underlying signal is stream of differentiated Diracs, in which case the conventional low-rank argument in |27| is 
not valid because it does not result in a Vandermonde decomposition structure of the associated Hankel matrix 
|. Therefore, Proposition |2. 1 1 will he used throughout the paper to justify the low-rankness. 

It is important to emphasize that Proposition |2.1| implies the following fundamental duality: 


T 


sparse signal low-ranked Hankel structured matrix, 











where T denotes the Fourier transform. Therefore, if some of the Hankel matrix elements are missing and one can 
measure samples only on the index set Q, the fundamental duality suggests a way to recover these elements, which 
is based on low rank matrix completion |[2^-|33|: 


(P) miUrngC" RANK,^(m) 

subject to Pn{in) = Fh(i 0 f) , 


(26) 


where 0 denotes the Hadamard product, and 1 and f denotes the vectors composed of discrete samples of i{uj) 
and f{u)), respectively. Note that by solving (P) we can obtain the missing data m{io) = l{u)f{(jj) in the Fourier 
domain. Then, the missing spectral data /(w) can be obtained by dividing by the weight, i.e. f{uj) = m{uj)/l{oj), 
when / 0. As for the signal /(w) at the spectral null of the filter f(w), the corresponding elements should be 
specifically obfained as sampled measuremenfs, which can be easily done in MR acquisition. 

Among various algorifhm fo solve mafrix completion problem (P), one of fhe mosf well-characferised approaches 


is a convex relaxation approach using fhe nuclear norm |291, p0| |, p2| , |331. More specifically, fhe missing k-space 
elemenfs can be found by solving fhe following nuclear norm minimizafion problem: 


(PI) rninmeC"! Il=^(m)||* 

subjecf fo Pn(m) = PQ(i 0 f) 


(27) 


where || • ||* denofes fhe mafrix nuclear norm. In Ibis case, we can further use fhe resulf by Chen and Chi |27| fo 
provide a performance guarantee. 


Theorem 2.2. [271 Suppose that the sampling index set is chosen randomly among [0,ni — 1] and |f2| = m. 
Then, if the number of k-space sample m is given by 

m > cipLiCsk\og^{ni) 


for some constant ci and an incoherence parameter pi, and 

Cs = max{ni/K, ni/(ni — k + 1)} 


(28) 


(29) 


then the perfect recovery of the missing spectral components is possible by solving (PI) with a probability exceeding 
1 - nj-^. 

Nofe fhaf fhe annihilafing fiber size k corresponds fo fhe “mafrix pencil paramefer” in fhe specfral compressed 
sensing approach by Chen and Chi |27|, who also showed fhaf fhe incoherence paramefer pi in pS]) grows fo one as 


fhe annihilafing filter size k increases p7| . Hence, fhe sampling rale in ( [28| ) is nearly opfimal since if is proporfional 
fo fhe unknown sparsify level up fo log^(-) factor. This suggesfs an imporfanf observation: by reformulaling fhe 
compressed sensing problem as a low rank Hankel sfrucfured mafrix complelion problem in fhe measuremenf 
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domain, no performance loss is expected. 

III. ALOHA FOR Accelerated MRI 

Inspired by the theoretical finding in the previous section, this section will explain two realizations of the low- 
rank matrix completion approaches that are useful for MR applications. The first one is a pyramidal decomposition 
algorithm that can be used for wavelet domain sparse signals, and the second one is a generalization for multichannel 
parallel MRI. 


A. Recovery of Wavelet Sparse Signals 

In fact, the signal / in Q can be sparsified using wavelet transform as demonstrated in |25|, | [26| . Since the 
sparsity in wavelet domain have been the main interest in the existing compressed sensing theory, we are particularly 
interested in this analysis approach. Interestingly, the resulting ALOHA framework can have very unique pyramidal 
structure that allows computational efficient and noise robust implementation of a low rank Hankel matrix completion 
algorithm. 


First, we assume that we have some real-valued “L-compatible” generalized wavelet [251, |26| which, at a given 
scale s, are defined as ipsix): 

ips{x) = L*(j)s{x) . (30) 


Here, L* is fhe adjoinf operafor of L and fs is a some smoothing kernel at s-scale. In dyadic wavelet decomposition. 


the s-scale wavelet is given by ipsix) = where 'ipo{x) is often called the “mother wavelet” | j^ . 

We further assume that the L-spline model in Q is cardinal, i.e. the knots positions {xj}^ZQ of the sparse driving 
signal ([T]) are on integer grid. 

Then, using (^i, the wavelet transform of the signal f{x) in ([^ at the s-scale is given by 


if, 'fsi- - x)) 


(/,L>,(.-x)) 

(L/, <Psi- - x)) 

-x)} = (f^*w)(x) 

k-1 

J2cjfs(x - Xj) 
j=0 


(31) 


where (/>g(x) = fs(—x) is the reversed version of fg, and the last equality comes from the definition of w in 0- 
For the case of “wavelet sparse” signals, their discrete wavelet transform (DWT) coefficients, which corresponds 
to the discrete samples of ( [3T| ) with the unit sampling interval, should be non-zero on sparse index set. Hence, if 
we choose maximally localised smoothing function fg, then ( [3T] ) can be made as sparse as possible. Ideally, if the 
(/>q(x) has the support size of 1, the sparsity of the 0-th scale DWT coefficients is equal to k, which achieves the 
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minimum sparsity level. In fact, this can be achieved using the triangle basis function p6| : 

0 < X < 1 

(/>o(x) =/3+(2x), where /3+(x) = 


X, 

2 — X, 1 < X < 2 

0, otherwise 


(32) 


which has the support size of 1; hence, it retains the sparsity level of the DWT coefficient to be same as k, which 
corresponds to the number of Diracs in the underlying sparse driving signal w{x). Furthermore, for the case of TV 


signals i.e. L = ^, the corresponding wavelet transform is indeed Haar wavelet, because we have 


1, 0 < X < ^ 


i/joix) = L*(j)o{x) = < -1, \<x <l 


(33) 


0, otherwise 


In general, the support size of (t)g{x) with respect to </>o(x) in (^i increases to 2®, so the maximum sparsity 
level at the s scale is up to 2®/c. To maintain the same sparsity level across the scale, we therefore use the dyadic 
wavelet decomposition, whose DWT coefficients are downsampled as 


f[l] ■■= 


(34) 


Due to the decimation by 2®, it is easy to see that the resulting sparsity level is now reduced to k. Then, the 
corresponding Fourier spectrum for the down-sampled signal is given by 

/®(e-):= 


(35) 

(36) 


where we use the scale dependent wavelet property 1351: 


and the righthand side terms of (361 come from spectral copies due to the downsampling. 


Because the number of nonzero coefficient /®[/] is at most k, then the form of the righthand side equation ( [35] ) 
is similar to ([T0|), so annihilating filters and the corresponding low-rank Hankel matrix can be found. However, 
one technical issue is that the righthand side of ( |3^ has aliased copies of k-space data weighted by wavelet 
weighting. Therefore, the direct application of low rank matrix completion may not work. Accordingly, to deal 
with 27r/2® repeating spectral structure, we propose a pyramidal decomposition of the Hankel structured matrix. 
More specifically, if the spectral component of /(w) is estimated for |a;| > 7r/2®, then we can constitute a residual 
signal recursively: 


fr{uj) : = 


0, |w| > 7r/2® 

otherwise 


(37) 
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(a) 


k -t Pyramidal steps 
(3 levels) 



2 scale ALOHA 

to! (10-') 

1 scale ALOHA 

tol (lO-'O 


low 


t high freq. 


scale ALOHA 

tol (10-0 


s scale ALOHA 



Coil 2 


■"n 


Coil 2 


(b) 

Fig. 1: ALOHA implementation using pyramidal decomposition. Construction of Hankel matrices from (a) kx-ky 
data by assuming that 2-D dydadic wavelet transform of images is sparse, and (b) k — t subsampled data by 
assuming that dynamic images can be sparsified using spatial wavelet and temporal Fourier transform. In the box, 
each reconstruction unit at the s-scale (the s-scale ALOHA) is illustrated, which consists of five steps: (1) the 
k-space region extraction, (2) k-space element by element weighting using s-scale wavelet spectrum 
(3) low-rank Hankel matrix completion, (4) k-space unweighting by dividing the interpolated k-space data using 
the s-scale wavelet spectrum, and (5) the k-space data replacement using the interpolated data. Note that for the 
case of dynamic MRI, one dimensional weigting is required along the phase encoding direction, whereas 2-D 
weighting is necessary for the case of static imaging. In the figure, ni denofes fhe ambient space dimension, and 
corresponds to the pseudo-inverse operation that takes the average value from the Hankel matrix and putting 
it back to the original k-space domain. The color coding in the Hankel structure matrix indicates the values of 
weighting. 


Then, for the corresponding discrete sample /*[(] := {fr^ipsi' — x))\^^ 2 h with fr = T ^{fr}, we have 

1 










( 38 ) 
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because its high frequency content will no more aliasing at a frequency inside |t<;| < 7r/2®. Then, due to the 
sparsity of there exists a filter that annihilates Consequently, the s-scale ALOHA problem can 

be formulated as 


min^gC"i/2° RANK,^(m) 
subject to Po(m) = Pq(L 0 f"" 


(39) 


where 


f[n] = f{uj) 


{jj=n'K ln\ 


and the weighting vector F is composed of following spectral samples 


f*[n] = F(a;) 


1 


(2*u;) 


(40) 


(41) 


u=n7Tlni 


LO=n7vlni 

where n = —ni/2*"’'^, • • • ,ni/2^'''^ — 1 and the superscript * denotes the complex conjugate. Note that sampling 
interval for 'ip*{co) widens by factor of 2 for each successive scale to construct the k-space weighting vector F. 
Consequently, we have a pyramidal decomposition as shown in Fig. [TJ and the ALOHA ( |^ should be solved 
from the lowest scale, i.e. s = 0, up to the highest scale (see caption of Fig. [T] for more discussion of the figures). 
Because sparsity can be imposed only on the wavelet coefficients, the low frequency k-space data that correspond 
to the scaling function coefficients should be acquired additionally during MR data acquisition. This information 
as well as the annihilating filter size then determines the the depth of the pyramidal decomposition as will be 
discussed in detail later. 

There are several advantages of using wavelet approaches compared to the direct operator weighting in ALOHA 
framework. First, the pyramidal decomposition structure for low rank matrix completion can significantly reduce 


the overall computational complexity. Moreover, it has been observed in |36| that the wavelet approach is more 
robust for noise and the model mismatch, which is also consistently observed in ALOHA framework as will be 
discussed later. 

Note that the proposed pyramidal scheme is an algorithm that is designed to approximately decouple the aliasing 
components in ( |3^ , and we do not claim the optimality of the proposed method. However, it is worth to mention 
that the optimality of a similar fine-to-coarse scale wavelet coefficients reconstruction was recently shown in Fourier 


compressed sensing problem |37|. 

B. Generalization to Parallel MRI 

Beside the annihilation property originating from sparsity in the transform domain, there exists an additional 
annihilation relationship that is unique in parallel MRI. The relationship we described here has similarity to SAKE 
and P-LORAKS when the image itself is sparse (even though these methods were inspired by the null space or 
GRAPPA type relationship), and our contribution is its generalization to transform domain sparse signals. 
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Specifically, in pMRI, the unknown image gi{x) from the i-th coil can he represented as 

gi{x) = Si{x)f{x), i = 

where Si{x) denotes the i-th coil sensitivity map, /(x) is an unknown image, and r denotes the number of coils. 
Then, for TV signals (i.e. L = ^), we have 

Lgi{x) = Si(x)Lf{x) + f{x)Lsi(x) , 

whose Fourier transform is given hy 

l{uj)gi{uj) = Si{uj) * {i{uj)f{u})) + f{uj) * (i{u})si(u)) . (42) 

Then, we can show the following inter-coil annihilating filter relationship: 

Lemma 3.1. Suppose that h{uj) is the annihilating filter for f{oj) such that {h * f){uj) = 0. Then, we have 

h{uj) * Sj{uj) * (J{uj)gi{uj)^ - h{uj) * sfiuj) * (j{u})gj{uj)^ =0, t 7 ^ j. (43) 

Proof. Using ( |42l ), we have 

h{ijj) * Sj{oj) * (j{uj)gi{Lo)^ = h{uj) * Sj{uj) * Si{u}) * {l{uj)f{u)) + h{Lo) * Sj{uj) * f{uj) * {l{u})si{Lo)) 

= h{uj) * Sj{uj) * Si{uj) * {l{uj)f{u!)) (44) 

where the second term vanishes thanks to {h* f){Lj) = 0. Similarly, we have 

h{uj) * sfiu}) * (l{u})gj{u})) = h{u}) * Si{uj) * Sj{uj) * {l{u})f{uj)) (45) 

□ 


Because Eqs. (44i and (451 are identical, Eq. (431 can he shown. Q.E.D. 


Therefore, to exploit the inter-coil annihilating property, Hankel matrix from each channel is stacked side hy side 
to form a matrix y: 


3^ = [jf'(iogi) ••• jf'(i0g^)] G, 


(46) 


Here, the augmented matrix y in (461 has the following rank condition: 


Proposition 3.2. For the concatenated Hankel matrix given in ( |46[ ), we have 

ranky < . 

where k and r denotes the sparsity and the number of coils, respectively. 


( 47 ) 


Proof. Since Eq. (431 holds for every pair among r-channels, it is easy to show that 


= 0 
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where 5i is defined recursively as follows: 

5.-1 ^ 




V. 

-v.-i 

Vl+l Vi+2 


VC 


-Vt 


-Vi 


5t+l 


-Vi 


where Vj denotes the discrete samples of h{ijj) * Si{uj). Accordingly, 


dimNUL(3^) > rank(5i) — ( 2 ) ~ ~ l)/^- 


Furthermore, Proposition 2.1 informs us that M’i} 0 gj) in y has rank at most k. Therefore, we have 

. , r(r-l) r(2k-r + l) 

rank 3^ < kr -= - . 


(48) 


(49) 


(50) 

□ 


Due to the low-rankness of the concatenated matrix, the multichannel version of the ALOHA can he formulated 
as 


min{m,}Ci rank [,9^(mi) • • • J^(m.)] (51) 

subject to Po(mi) = Pn(i 0 g*), z = 1, • • • , r . 


In addition, it is worth to emphasize the meaning of the rank condition of the concatenated matrix. Because the 
rank of y is at most r{2k — r + l)/2 and there exist additional degrees of freedom that belong to the amplitudes 
of Diracs in the transform domain, the total degrees of the freedom for the FRI signal are at most r{2k — r + 1) 
In parallel MRI, if m denotes the number of k-space sampling locations, then k-space data are sampled 
simultaneously from r-coils and the total number of k-space samples becomes mr. Since the number of the samples 
should be larger than the degree of the freedom, we have 


mr > r{2k 


r + l) 


k < 


m + r — 1 


(52) 


This result has very important geometric meaning. Suppose that we attempt to address parallel imaging by exploiting 
the joint sparsity. Using the notation in ([1]), this could be addressed by the following multiple measurement vector 
(MMV) problem 


min ||P||o 
F 

subject to G = AF (53) 


where G = [gi, • • • , g.] and P = [fi, • • • , f.] with fj = 
number of non-zero rows. Then, it was shown in our 


= [S'i]f for a given sensitivity map [Si], and || • ||o denotes the 
previous work |22| and others | [2^ , p4| that if P* satisfies 
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AF^ = G and 


\F4o < 


spark(yl) + RANK(G) — 1 


(54) 


then F* is the unique solution for G = AF, where spark(^) denote the smallest number of linearly dependent 
columns of A. For randomly chosen Fourier samples, we have spark(^) > m + 1. Moreover, RANk(G) = r. 

if the unknown sparsity level is k, i.e. ||F*||o = k. 


Hence, ( [54| ) is equivalent to 

Therefore, similar to the fact that single coil ALOHA does not lose any theoretical optimality compared to 
the single coil compressed sensing approaches, our ALOHA formulation for pMRI using ( [ST] ) may fully exploit 
the multi-channel diversity from parallel acquisition. Moreover, because the k-space interpolation formulation of 
ALOHA is derived based on the information of the singularity of the signals, the resulting reconstruction results 
for both single and multi-channel formulation exhibit superior performance as will be discussed in experimental 
sections. 

IV. Implementation Details 
A. 2D Hankel Structured Matrix Construction 

In 3D imaging or dynamic acquisition of MR data, the readout direction is usually fully sampled and the other 
two encoding directions are under sampled. Thus, this section presents an explicit way of constructing a 2D Hankel 
structured matrix. Specifically, if h[n, m] is a pi x qi size 2D annihilating filter, then the corresponding annihilating 
filter relation is given by 

Pi-i (ji-i 

{h* f)[n,m] = (55) 

i=0 j=0 

for all n,m £ Q. Let ni x mi k-space data matrix be defined by 

/[0,0] ••• /[0 ,mi-l] 

F := : : 

_/[ni-l,0] ••• /[ni-l,mi-l] 

= [fo • • • fmi-l] 

Similarly, we define pi x qi annihilating filter matrix H. Then, by removing the boundary effect from the 2D 
convolution as shown in Fig. [^a), the 2D annihilation property should hold only inside of the domain and (551 
can be equivalently represented as 


A^(F)I) = 0, 

where a 2-D Hankel structured matrix Jif{F) is constructed as 


(56) 


^(fo) 

J^(fi) 






Ji-iiqj 


( 57 ) 
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(d) 

Fig. 2: (a) Area where annihilation property holds. Various ways of constructing block Hankel matrices: (h) ALOHA, 
(c) SAKE, and (d) LORAKS. In (d). Nr denotes the number of neighborhood pixels. 


with G Pi+i)xpi given by 


/jO,j] 


/[2,j] 


f[ni-pi,j] f[ni-pi + l,j] 

and the annihilating filter vector is given by 


f\pi - 1,2] 

flm -i,j] 


l) = Yec{H) 


(58) 


where the overline denotes an operator that reserves the order of a vector. Using this, we can construct an augmented 
matrix y in ( |46l ) from r-channels. 

The augmented matrix structure Jif{F) illustrated in Fig. |^b) is similar to those of SAKE and EORAKS/P- 
LORAKS in Fig. |^c) and (d), respectively, with the following differences. Compared to SAKE, AEOHA stacks the 
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multi-coil Hankel matrices side by side. Unlike the SAKE and ALOHA, LORAKS uses i? = 4 {R: neighborhood 
size) neighbors according to the software manual provided by the original authors. Moreover, the most important 
novelty of the proposed method is the k-space weighting to construct the Hankel structured matrix. 


B. Hankel structured matrix completion algorithm 


In order to solve Eqs. (26) and (511, we employ an SVD-free structured rank minimization algorithm |38| with 


an initialization using the low-rank factorization model (LMaFit) algorithm |39|. This algorithm does not use the 


singular value decomposition (SVD), so the computational complexity can be significantly reduced. Specifically, 
fhe algorithm is based on the following observation 


= mm 
U,V:A=UV» 


2 -e 
F + 


(59) 


Hence, (26 1 can be reformulated as the nuclear norm minimization problem under the matrix factorization constraint: 


mm 

U,V:jr{m)=UV" 


2 -e 
F + 


subject to Eb(m) = Pn(f). 


(60) 


By combining the two constraints, we have the following cost function for an alternating direction method of 


multiplier (ADMM) step |41|: 


L([/,U,m,A) := ,{m) + U\\U\\l + \\V\\l) 


+|||J^(m)-C/U^ + A||| 


(61) 


where i(m) denotes an indicator function: 


r(m) = |°’ ifPo(m) = Pn(f) 

\ oo, otherwise 

One of the advantages of the ADMM formulation is that each subproblem is simply obtained from ( [6T] ). More 
specifically, and can be obtained, respectively, by applying the following optimization 

problems sequentially: 


minmi(m) + ^||Af'(m) - + A(”)||| 

1 iirni2 ^ M||^(m(^+i)) _ 

1 


(62) 


mmy 5II»/ ill. + f ||^(m("+i)) - + A^ ||2, 

and the Lagrangian update is given by 

A(^+1) = y{n+l) _ jj{n+l)y[n+l)H _|_ ^(n) 


It is easy to show that the first step in ( [62] ) can be reduced to 

,(’^+1) = ^u{n)y{n)H _ ^(n) | ^ 




(63) 
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where Pqo is a projection mapping on the set and corresponds to the Penrose-Moore pseudo-inverse mapping 
from our block Hankel structure to a vector. Hence, the role of the pseudo-inverse is taking the average value and 
putting it hack to the original coordinate. Next, the suhprohlem for U and V can he easily calculated hy taking the 
derivative with respect to each matrix, and we have 

, -1 


jjin+l) ^ ^ ( 3 ;(n-rl) ^(n)^ y(n) 

y{n+l) _ ^ ^^(n-l-1) _|_ jjiri+l) _|_ ^u{n+l)H 


(64) 


Note that the computational complexity of our ADMM algorithm is dependent on the matrix inversion in ( |64| ), 
whose complexity is determined hy the estimated rank of the Hankel matrix. Therefore, even though the Hankel 
matrix has large size, the estimated rank is much smaller, which significantly reduces overall complexity. 

Now, for faster convergence, the remaining issue is how to initialize U and V. For this, we employ an algorithm 


called the low-rank factorization model (LMaFit) |39|. More specifically, for a low-rank mafrix Z, LMaFif solves 
fhe following optimization problem: 


-||f71/^ — Z\\‘^p subject to Pi{Z) = P/(J^(f)) 


(65) 


and Z is initialized with and the index set I denotes the positions where the elements of JZ{i) are known. 

LMaFit solves a linear equation with respect to U and V to find their updates and relaxes the updates by taking 
the average between the previous iteration and the current iteration. Moreover, the rank estimation can be done 
automatically. LMaFit uses QR factorization instead of SVD, so it is also computationally efficient. Even though 
the problem ([65|) is non-convex due to the multiplication of U and V, the convergence of LMaFit to a stationary 


point was analyzed in detail |39|. However, the LMaFit alone cannot recover the block Hankel structure, which is 
the reason we use an ADMM step afterward to impose the structure. 


C. Reconstruction Flow 

As shown in Fig. [TJ the ALOHA framework is comprised with several major steps: pyramidal decomposition, 
k-space weighting, Hankel matrix formation, rank estimation, SVD-free low rank matrix completion, and k-space 
unweighting. Here, we will explain these in more detail. 

The pyramidal decomposition is performed as follows. First, in static MR data acquisition illustrated in Fig. [TJa), 
the kx — ky corresponds to the two phase encoding directions that are downsampled. Thus, the Hankel matrix is 
constructed from kx — ky data. After a k-space interpolation from a finer scale, the data at the current scale is 
defined to contain one-fourth of data around zero frequency from that of the previous scale. Second, in the case 
of dynamic MR imaging shown in Fig. [^b), kx samples from the readout direction are fully acquired, whereas 
the ky directional phase encoding are downsampled along the temporal direction t. Therefore, the data in ky — t 
space (or simply, k — t space) are downsampled, from which we construct a Hankel structure matrix. In pyramidal 
decomposition, after a k-space interpolation from a finer scale, fhe ky — t dafa in the current scale contains a half 
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of the data from that of the previous scale. Note that the wavelet decomposition is performed only along the spatial 
domain, so the pyramidal decomposition is only performed along ky direction. This construction of Hankel matrix 
is due to the ohservation that the dynamic signal is sparse in spatial wavelet and temporal Fourier transform domain 
m- See more details in our recent work p2| . 

In both cases, the estimated k-space data at the lower scale are used to initialize the low rank matrix completion 
algorithm at the current scale. This accelerates the convergence speed. Moreover, due to the additional chance 
of refining the estimates, more important k-space samples at the low frequency regions are refined furfhermore 
compared fo fhe high frequency k-space samples. Consequently, the overall computational burden of the low rank 
matrix competition algorithm is significantly reduced while the overall quality is still maintained. 

The k-space weighting is performed using wavelets. Specifically, based on our discussions on maximally localized 
smoofhing funcfion, we use a Haar wavelef expansion whose specfrum is given by 

= f (!^) exp(-i|) , (66) 

The corresponding k-space weighting af fhe s-scale is given by 

^ 1 /sin2*a;/4 

^ V 2*u;/4 



For fhe case of sfafic MRI in Fig. [TJa), we use 2-D weighting by assuming fhaf fhe image is sparse in 2-D dyadic 
wavelef fransform domain. Care needs fo be taken when applying the weighting to 2D Fourier domain because 
there are two frequency variables {ujx,i 0 y)- One could use a separable weighting l{ujx,i 0 y) = l{uJx)K^y)’ however, 
the resulting problem is that the missing k-space components along the frequency axis iOx = 0 or ujy = 0 cannot 
be recovered. Consequently, we applied the weighting sequentially along each axis, i.e. we solve ( |26l ) by applying 
l{iOx) first, which is followed by solving (|2^ with l{iOy). However, simultaneous weighting would be possible as 


demonstrated in a recent work for off-the-grid recovery of piecewise constant image |43|. For the case of dynamic 
imaging in Fig. [^b), one dimensional weighting along the phase encoding direction was applied as explained in 
detail in p2| . Finally, after the k-space interpolation, the k-space unweighting is done in k-space pixel-by-pixel by 
dividing the reconstructed value with ( |6^ . Note that (|6^ has zero value at the DC frequency. However, because 
we acquire the DC value as well as some of the low frequency k-space data, the problem of dividing by zero never 
happened. 

Because LORAKS is similar to the ALOHA without weighting, we use LORAKS as a reference to contrast why 
the proposed ALOHA framework has many advantages. The implementation of LORAKS was based on the source 
code available in the original author’s homepage, which requires manual setting of estimated ranks. We chose the 
rank for LORAKS that gave the best reconstruction quality. In Discussion, we also provided reconstruction results 
by ALOHA without weighting to control the confounds and confirm fhe imporfance of fhe k-space weighfing in 
consfructing Hankel mafrix. 
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We used TITAN GTX graphic card for graphic processor unit (GPU) and i7-4770k CPU and the codes were 
written in MATLAB 2015a (Mathwork, Natick). To accelerate the algorithm, most part of the MATLAB codes 
were implemented using Compute Unified Device Architecture (CUDA) for GPU. 

D. MR Acquisition and Reconstruction Parameters 

To assess the performance of ALOHA for single coil compressed sensing imaging, k-space raw data from an 
MR headscan was obtained with Siemens Tim Trio 3T scanner using balanced steady-state free precession (bSSFP) 
sequence. The acquisition parameters were as follows: TR/TE = 10.68/5.34ms, 208 x 256 acquisition matrix 
(partial Fourier factor 7/8, oversampling factor 50%), and number of slices is 104 with 2mm slice thickness. The 
field-of-view (FOV) was 178 x 220mm^. We used the central coronal slice. 

A retrospective down-sampling mask was generated according to a two dimensional Gaussian distribution and the 
data at the central 7x7 region around zero frequency were obtained additionally. This is equivalent to assume a 3D 
imaging scenario where the readout direction is fully sampled and the downsampling is done in the remaining 2-D 
phase encoding direction. Downsampling factors of four was used to generate sampling masks. However, data was 
obtained as partial Fourier measurements, so effective downsampling ratio was 4.13. The 2-D k-space weighting 
using ( [^ was used. The AFOHA reconstruction was conducted using the following parameters: three levels of 
pyramidal decomposition, and decreasing FMaFit tolerance values (5 x 10“^, 5 x 10“^, 5 x 10“^) at each level of 
the pyramid. In addition, an initial rank estimate for FMaFit started with one and was refined aufomafically in an 
increasing sequence, fhe annihilating filler size was 23x23, and Ihe ADMM paramefer was /r = 10^. 

For fhe compressed sensing approach, we used Iwo approaches wifh fhe same dala and fhe same sampling masks: 
(1) fhe sparsity in wavelet domain (which we denote (i-wavelet) 0, and (2) the split Bregman method for the total 


variation 144|. In case of /i-wavelet approach, we implemented an ADMM algorithm using wavelet domain sparsity. 


In addition, for comparison with the existing state-of-the art approach using Hankel structured matrix completion 
algorithm, FORAKS was compared. In particular, C-based FORAKS was chosen because it exploits the image 
domain sparsity. The parameters for the fi-wavelet and TV approaches were optimized to have the best performance 
in terms of the normalized mean square error (NMSE), which is defined by NMSE(x) = ||x — ylH/llylli’ where x 
and y denote the reconstructed and the ground-truth images, respectively. The FORAKS parameters were chosen 
manually to give the best reconstruction quality. 

To evaluate the performance of AFOHA in static parallel imaging, k-space raw data from an MR headscan 
was obtained with Siemens Verio 3T scanner using 2D SE sequence. The acquisition parameters were as follows: 
TR/TE = 4000/100ms, 256 x 256 acquisition matrix, and six z-slices with 5mm slice thickness. The field-of-view 
(FOV) was 240 x 240mm^, and fhe number of coils was four. Refrospecfively undersampled 2-D k-space dala al fhe 
acceleralion factor of eight were obtained according to a two dimensional Gaussian distribution in addition to the 
7x7 central region around zero frequency. The data from 4 receiver coils were used. For comparison, we used the 
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identical data and sampling masks for SAKE | [T6| , and SAKE with ESPIRiT |451. Note that GRAPPA | fT3| requires 
ACS lines, so with the additional 50 samples along ACS, the effective downsampling ratio was 4.785, which is not 
good as the eight time acceleration in other algorithms. SAKE and SAKE with ESPIRiT are both low rank matrix 
completion algorithms for Hankel structured matrix collected from the whole k-space data. However, SAKE with 
ESPIRiT reduces the computational burden of the original SAKE by performing low rank matrix completion only 
for the 65 x 65 central region with 7x7 filter, after which coil sensitivities are estimated using the reconstruction 
data. The estimated coil sensitivities are used to estimate the remaining k-space missing data through ESPIRiT | |45| . 
The parameters for SAKE and SAKE with ESPIRiT were chosen such that they provided the best reconstruction 
results. The parameters for the AEOHA are as follows: three levels of pyramidal decomposition with decreasing 
EMaEit tolerances (10“^, 10“^, 10“^), and 7x7 annihilating filter. The same EMaEit rank estimation strategy and 
ADMM parameter used for single coil experiments were employed. We generated the square root of sum of squares 
(SSoS) image from multi-coil reconstructions. 

We also validated the performance of AEOHA for accelerated dynamic cardiac data in the k — t domain. A cardiac 
cine data set was acquired using a 3T whole-body MRI scanner (Siemens; Tim Trio) equipped with a 32-element 
cardiac coil array. The acquisition sequence was bSSEP and prospective cardiac gating was used. The imaging 
parameters were as follows: EOV= 300 x 300mm^, acquisition matrix size= 128 x 128, TE/TR = 1.37/2.7ms, 
receiver bandwidth = 1184Hz/pixel, and flip angle = 40°. The number of cardiac phases was 23 and the temporal 
resolution was 43.2ms. The k-t space samples including four lines around zero frequency were retrospectively 
obtained at the reduction factor of eight according to a Gaussian distribution. Eor comparison, k-t EOCUSS Q, 
EORAKS p7| , and SAKE was used. In case of EORAKS (single coil) and SAKE (multi coil), these algorithms 
were applied to k-t domain for dynamic reconstruction. The parameters in k-t EOCUSS, EORAKS, and SAKE, were 
selected to give the best NMSE values. Eor the AEOHA in single coil data, the following parameters were used: 
three level of pyramidal decomposition only along the phase encoding direction, decreasing EMaEit tolerances 
(10“^, 10“^, 10“^) at each scale, and 17x5 annihilating filter. The same EMaEit rank estimation strategy and 
ADMM parameter was // = 10. The k-space weighting in Eq. ( [66l ) was applied only along the phase encoding 
direction. 

Next, we investigated the synergetic improvement of dynamic imaging from multi-channel acquisition. Eour 
representative coils out of 32 were used. The reason we chose only four coils was to verify that the generalised 
AEOHA can maximally exploit the multi-channel diversity even with the small number of coils. The four coils 
were chosen such that it covers every area of images. The same four coils were used for all algorithms for fair 
comparison. In the AEOHA, the annihilating filter size was 7x7. The same EMaEit rank estimation strategy and 
ADMM parameter used before were employed. After the reconstructions of k-space samples, the inverse Eourier 
transform was applied, and the SSoS images were obtained by combining the reconstructed images. 
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Wavelet (x4) TV (x4) LORAKS (x4) Proposed (x4) 

ORIG NMSE 1.912e-02 NMSE 1.581e-02 NMSE 1.702e-02 NMSE 1.349e-02 



Fig. 3: Comparison with -wavelet compressed sensing, TV compressed sensing, LORAKS, and the proposed 
method at four fold acceleration factors. The data was acquired from a single channel coil. The first row shows 
reconstructed images, and the second row shows difference images between the ground-truth and the reconstructions, 
and the third row shows the magnified views of distorted regions in the reconstructed images. The last row shows 
the difference images in the magnified views. 


V. Results 

A. Static MR experiments 

Reconstrucfed results from a single coil hrain data are shown in Fig. with the NMSE values. From the NMSE 
values, we observed that the performance AEOHA was quantitatively superior to the performance of Zi-wavelet 
and TV based compressed sensing approach. The reconstruction results by AEOHA has less perceivable distortion 
compared to those of Zi-wavelet and TV approaches. This can be easily observed from the difference images in 
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GRAPPA (x8) SAKE (x8) SAKE + ESPIRiT (x8) Proposed (x8) 


ORIG NMSE 1.242e-02 NMSE 2.843e-03 NMSE 2.644e-03 NMSE 2.210e-03 



Fig. 4: Parallel imaging results using GRAPPA, SAKE, SAKE with ESPIRiT and the proposed method at 8 
fold acceleration. The second row shows the difference images. Areas with systematic artefacts are indicated hy 
yellow arrows. Note that GRAPPA requires ACS lines, so with the additional 50 samples along ACS, the effective 
downsampling ratio was 4.785. 


the second and the fourth rows of Fig. In the case of ii-wavelets and TV, structural distortion around the image 
edges was easily recognizable. In the last row of Fig. the edges were reconstructed accurately hy ALOHA. On 
the other hand, the contrast between grey matters and white matters in -wavelets and TV reconstruction were 
significantly distorted compared with that of ALOHA. The LORAKS reconstruction appeared better than that of 

-wavelets and TV reconstruction, but there were still remaining errors around the edges and the NMSE value was 
much higher than that of ALOHA. 

Next we compared our parallel imaging results with those of the existing approaches for additional multichannel 
brain data set. The NMSE results in Fig. showed that ALOHA was most accurate. From the difference images 
at the second row of Fig. we observed that proposed method provided reconstruction results more accurately 
than other algorithms. In SAKE, the structures were distorted around the inner skull and the boundaries of tissues. 
In SAKE with ESPIRiT, overall reconstruction errors were higher than those from ALOHA and there were still 
remaining errors around the skull. The reconstruction time was 22.2sec with our preliminary GPU implementation 
of ALOHA, which attained a speed-up factor of 5 compared to CPU implementation. On the other hand, the 
computational time for MATLAB version of GRAPPA, SAKE, and SAKE-i-ESPIRiT were 9s, 320.4s, and 21.6s, 
respectively. 

B. Dynamic MR experiments 

Using the sub-sampled k-space data at the acceleration factor of eight, the average NMSE values of k-t EOCUSS, 
EORAKS, and AEOHA, were 1.616 x 10“^, 1.363 x 10“^, and, 1.151 x 10“^, respectively. The sub-sampled data 
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(a) Single coil reconstruction results 


ORIG 



k-t FOCUSS (x8) 
NMSE 1.616e-02 


LORAKS (x8) 
NMSE 1.363e-02 


Proposed (x8) 
NMSE 1.151e-02 
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(b) Parallel imaging results from 4 coils 

Fig. 5: Reconstruction results from 8 fold accelerated k-space data using (a) single coil and (b) four coils data 
set. Purple lines denote the regions corresponding to y-t cross sections that are magnified along temporal axis. The 
second rows in both (a) and (b) show the difference images between the ground-truth and the reconstructions. 



k-t FOCUSS (x8) 
NMSE 8.750e-03 


SAKE (x8) 
NMSE 4.983e-03 


Proposed (x8) 


was collected according to a Gaussian distribution and included the four center lines around zero frequency. The 
average NMSE values were calculated using all temporal frames. These results confirmed fhat fhe proposed method 
outperformed k-t FOCUSS and LORAKS. As shown in Fig. |^a), the temporal profile (indicated as a broken purple 
line) of the proposed reconstruction provided more accurate structures especially in the systole phase that were 
comparable to the true one, whereas the temporal variation in the k-t FOCUSS and LORAKS reconstruction became 
smoother and more blurry along the temporal dimension. 

The NMSE values of the parallel dynamic imaging results from k-t FOCUSS, SAKE, and the proposed method 
using four coil k-space data were 8.75 x 10“^, 4.983 x 10“^ and 3.436 x 10“^, respectively, which quantitatively 
showed that the proposed method outperformed k-t FOCUSS and SAKE. Reduced residual artifacts were perceivable 
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in the ALOHA difference images in Fig. |^b). Moreover, the temporal profiles of the proposed reconstruction 
showed more accurate structures which were comparable to the true one, whereas the dynamic slice profile from 
k-f FOCUSS and SAKE showed smoofher and more blurry fransifion. In parficular, fhe proposed algorifhm resulfed 
in more accurate reconstructions of dynamic changes at the heart wall in the systole phase as shown in Fig. |^b). 

VI. Discussion 

A. Effects of k-space weighting 

Figures [^a)(b) illustrate the effects of wavelet weighting schemes. In order to also demonstrate the sensitivity 
with respect to annihilating filter size, we calculated the NMSE values by changing the filter size. Eurthermore, 
to decouple the confounds originated from different implementation of C-based EORAKS and the AEOHA, the 
experimental results were generated using the same AEOHA framework with identical parameter setting, except for 
those related to the weighting. Eigure |^a) showed the single channel brain reconstruction results at the acceleration 
factor of 6. With weighting, the AEOHA reconstruction was conducted using the following parameters: three levels 
of pyramidal decomposition, and decreasing EMaEit tolerance values (5 x 10“^, 5 x 10“^, 5 x 10“^) at each level 
of the pyramid. In addition, an initial rank estimate for EMaEit started with one and was refined automafically 
in an increasing sequence, and fhe ADMM paramefer was p = 10^. Nofe fhaf fhese paramefers were same wifh 
fhose for Q. Eor fhe case of non-weighfed implementation of AEOHA, the other parameters are exactly the same 
except the weighting. The results showed that the NMSE values of the weighted AEOHA are consistently better 
than those of unweighted AEOHA regardless of the annihilating filter size. Moreover, the NMSE values were not 
sensitive to the annihilating filter size for the cases of the proposed AEOHA, whose NMSE values converged. The 
reconstruction results from two implementation at the minimum NMSE values (marked as stars in Eigures [^a)(b) 
) were illustrated, which again clearly showed that the residual errors of the weighted AEOHA is significantly 
smaller than the unweighted version of AEOHA. 

Similar results were obtained from the dynamic cardiac imaging data in Eigure [^b). With weighting, the 
AEOHA reconstruction was conducted using the following parameters: three levels of pyramidal decomposition, 
and decreasing EMaEit tolerance values (10“^, 10“^, 10“^) at each level of the pyramid. In addition, an initial rank 
estimate for EMaEit started with one and was refined automafically in an increasing sequence, and the ADMM 
parameter was p = 10. The reconstruction parameters for the unweighted implementation of AEOHA were exactly 
the same except the weighting. The results showed that the unweighted AEOHA is very sensitive to the annihilating 
filter size, which showed the divergent behavior as the filter size increases. However, the proposed weighted AEOHA 
exhibited the convergent behaviors. This result clearly confirmed Proposition |2T] saying that the low rank structure is 
invariant as long as the annihilating filter size is bigger than the transform domain sparsity level. The reconstruction 
results at the minimum NMSE values clearly showed that the proposed approach provided more clearly transition 
between diastole and systole phases. 
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(a) single coil brain data (k^-ky) x6 acceleration 



Without weighting With weighting 



size of annihilating filter 



(b) single coil cardiac data (k-t) xlO acceleration 

Fig. 6: The effect of weighting and annihilating filter size in ALOHA implementation, (a) Single coil hrain data 
results, and (h) single coil dynamic cardiac imaging results. The results clearly showed that the weighting is 
necessary for ALOHA to exploit the transform domain sparsity. On the other hand, the unweighted implementation 
exhibited divergent behavior with increasing filter size. 


B. Hyper-Parameter Estimation 

In order to show the importance of the pyramidal decomposition. Fig. [ 7 ] plots the computational time versus 
NMSE values of single coil brain data by changing the number of decomposition levels. As discussed before, the 
maximum decomposition level is determined by the acquired low frequency components and the annihilating filter 
size. In this data set, k-space dimension was 208 x 512 and the annihilating filter size was 23 x 23, and 7x7 k-space 


data around zero frequency were acquired. By converting the Hankel matrix dimension in ( pT] ) for the s-scale, the 
number of rows should not be smaller than the number of columns; so we should impose constraint 


ni/2^ - Pi + 1 = 208/2^ - 23 + 1 > 23, 
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Fig. 7: Reconstruction NMSE values with respect to different decomposition levels. The data is from single channel 
static MR reconstruction results in Fig. 


where s is the scale. This provides s < 2, and the maximum scale becomes 2. Fig. showed that the performance 
gain increases as a scale increases; as expected, for s > 3, the performance improvement was negligible. Therefore, 
the maximum scale was determined as s = 2. 

The other important reconstruction parameters include the size of the annihilating filter, the number of iterations, 
and tolerances used in LMaFit algorithm. Recall that the annihilating filter size corresponds to the matrix pencil 


size in sensor array signal processing |34|, and it should be set larger than the sparsity level of the transform 
coefficients. In fact, by considering the expression of c* in ( |29l ), the optimal annihilating filter size for 2D data 
is ni/2 X mil2, where ni x mi denotes the full k-space data size. Flowever, such large annihilating filter size 
introduces significant computational burden, so we tried to reduce the filter size as long as the image quality is 
not degraded. Based on extensive experiments, we found that in single coil image, the filter size should be set 
larger than that of parallel imaging, because the annihilating filter size is solely determined by the sparsity level. In 
parallel imaging, there exist additional annihilating filters from the intercoil relationship, so the annihilating filter 
size for each Hankel matrix construction can be set smaller than that of a single coil imaging. 

Finally, the tolerance level for LMaFit, which corresponds to the fitting accuracy, plays key role in determining 
the initial rank estimate. The initial rank estimate need not be close to the exact rank, but it was used to define fhe 
dimension of U and V in ADMM. We found fhaf fhe folerance level could be defermined by considering differenl 
noise contribufions in k-space dafa. Specifically, higher frequency components are usually contaminated by higher 
level of noises compared to the low frequency k-space data, so the LMaFit fitting accuracy need not be enforced 
strictly. This was the case when LMaFit was applied at a lower scale, since high frequency k-space data are more 
weighted and noises were boosted. On the other hand, more accurate fitting is required for higher scale data where 
the lower frequency k-space data are more weighted. Consequently, we chose decreasing values of tolerances per 
scale for in-vivo experiments. 
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ALOHA ALOHA ver. 

ORIG NMSE 2.210e-03 NMSE L428e-02 



Fig. 8: Dependency of concatenation direction of Hankel matrix for parallel MRI: (middle) side by side augmentation 
of Hankel matrices, and (right) vertical augmentation of Hankel matrices. Side-by-side augumentation provided the 
superior reconstruction results. 


C. Concatenation Direction of Hankel Matrices for Parallel MRI 

When generalizing the Hankel matrix model from single to multiple channels, we chose to stack the Hankel 
matrices corresponding to individual channels side by side. One could suggest stacking them one on top of another, 
which would represent a different matricization. Based on our discussion in Section]^ in order to stack the Hankel 
matrix on top of another, the Hankel matrices should share the same annihilating filter. However, for the TV 


signal, the structure in (^i showed that it is not possible to find a common annihilafing filter hciijj) such that 
hc{oj) * H{uj)gi{uj) j = 0, Vi. On the other hand, the relationship in Eq. ( |4^ clearly showed that there are inter-coil 


annihilating filters that can cancel out each combinations of the k-space signals. Therefore, to utilize ( |43] ), we should 
stack the Hankel matrix size by side. 

To confirm fhe theoretical findings, we compared the reconstruction results from the two different stacking of 
Hankel matrices. As shown in Fig. the vertical augmentation of Hankel matrices provided inferior reconstruction 
results compred to the side by side augmentation of Hankel matrices. This again confirmed fhe importance of the 
concatenation direction and the analysis by Proposition |3.2| 


D. Further Extensions 

Note that this work is a first step towards unifying sparse and low rank models, so there are many rooms for 
improvement. In particular, it could be extended for other measurement and sparsity models. For example, while 
Eq. ( |2^ is based on a data equality constraint that focuses on standard Cartesian image reconstruction, we could 
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convert (P) in ( [26l ) to deal with the cases of general non-Cartesian imaging and/or noisy measurements: 

(P') minmgC" RANK^(m) (67) 

subject to ||Pn(m) - AQ(i©f)|| < (5 , 

where m denotes the spectrum data on cartesian grid, Aq, denotes a linear mapping that interpolates non-cartesian 
data iQf to the nearest cartesian grid index and the noisy level S is determined hy the gridding or measurement 
noises. The minimization problem (P') can be also addressed using an SVD-free ADMM approach similar to ( |6T] ), 
except that the indicator function-based data fidelity term is changed to I 2 data fidelity term. The detailed algorithm 


is described in |28|. 


In addition, ALOHA could be extended for reconstruction models that incorporate various system non-idealities. 


For example, in our recent work |46|, we showed that EPI ghost artifacts that originate from off-resonance related 
inconsistencies between odd and even echoes can be removed by exploiting that the differential k-space data 
between the even and odd echoes is a Fourier transform of an underlying sparse image. Specifically, we can 
consfrucf a rank-deficienl concafenafed Hankel sfrucfured mafrix from even and odd k-space dafa, whose missing 
dafa can be inferpolafed using ALOHA. However, fhe exfension of ALOHA for general off-resonance correcfed 
MR reconstruction is still an open problem, which needs a further investigation in the future. 

While this paper assumes that signals can be sparsified using dyadic wavelet ransforms, in general, signals can 
be more easily sparsified using non-decimafed redundanf wavelefs. In fhis case, fhe corresponding non-decimafed 
discrete sample af fhe s-scale is given by /^[/] := — x))\^^i whose specfrum can be represented by 

1 


I 


—iul 




V^*(2®a;)/(w). 


( 68 ) 


Because fhere exisfs no aliasing componenfs due fo fhe lack of downsampling, fhe construction of weighted Hankel 
matrix is much simpler than the dyadic wavelet transform. However, one of the potential downsides is that the 
number of non-zero wavelet transform coefficients /*[/] increases up to 2^k. Therefore, more k-space measurements 
are required to recover k-space data corresponding to the coarse scale wavelet coefficients. One potential solution 
would be to utilize the recovered high frequency k-space samples as additional measurements for interpolating the 
coarser level k-space samples. Note that this is different from dyadic wavelet case which discards higher frequency 
k-space samples in recovering low frequency k-space samples. However, the efficacy of fhis proposal needs fo be 
evaluafed sysfemalically, which is beyond fhe scope of fhis paper. 

Ofher fhan wavelef or TV represenfafion, modern compressed sensing approach deals wifh advanced sparsifying 
Iransforms such as pafch-based mefhods, dicfionary learning, Markov penalfies, and so on. The exfension of ALOHA 
for such model would be imporfanf and rewarding, which requires more extensive invesfigafion in subsequenf 
research. Recenfly, ALOHA was successfully used for MR paramefer mapping | [42| . Considering recenf success 
of direcf paramefer mappings |47|, | |48| , fhe exfension of ALOHA for direcf paramefer mapping would be very 
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interesting, which deservers further investigations in the future. 

VII. Conclusion 

In this paper, we proposed a general framework for annihilating filter based low-rank Hankel matrix approach 
(ALOHA) for static and dynamic MRI inspired by recent calibration-free k-space methods such as SAKE and 
LORAKS/P-LORAKS. Because natural images can be much more effectively sparsified in the transform domains, 
we generalized the idea to include signals that can be sparsfied in the transform domains. Our analysis showed 
that the transform domain sparsity can be equivalently represented as low-rank Hankel structured matrix in the 
weighted k-space domain, whose weighting function is determined solely by the transform, not by the data. In 
particular, when signals are effectively sparsified in dyadic wavelet transform, the corresponding low rank Hankel 
matrix completion problem can be implemented using a pyramidal decomposition, which significantly reduces the 
overall computational complexity and improves the noise robustness. For parallel imaging data, we verified that by 
stacking Hankel matrix from each coil side by side, we may fully exploit the coil sensitivity diversity thanks to 
existence of inter-coil annihilating filters. 

Reconstruction results from single coil static MR imaging confirmed that the proposed method outperformed the 
existing compressed sensing framework with TV regularization. We further demonstrated superior performance of 
the proposed method in static parallel MR imaging even without calibration data. Furthermore, the algorithm was 
successfully extended to dynamic accelerated MRI along k-t domain with both single coil and multi coil dynamic 
MR data. Therefore, we concluded that the proposed algorithm was very effective in unifying the compressed 
sensing and parallel MRI. 
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